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ABSTRACT 

The  mathematical  tools  for  analyzing  the  two-dimensional  and  axlsymmetric  flow  fields 
of  a  supersonic  stream  in  which  chemical  reactions  are  occurring  have  been  available 
for  several  years.  The  numerical  computations  of  these  flow  fields  were  not  considered 
practical  until  the  introduction  of  high-speed  computers  and  the  solution  cf  such  a  problem 
for  the  one-dimensional  case  had  been  completed.  Recently,  a  computer  program  was 
developed  that  applied  the  numerical  solution  of  the  method  of  characteristics  with  finite 
rate  chemical  reactions.  This  program  enables  one  to  solve  two-dimensional  and  axisym- 
metrlc  reactive  flow  fields.  This  report  discusses  the  difficulties  encountered  and  their 
solutions  in  the  development  of  the  computer  program  for  the  two-dimensional  and  axisym- 
metrlc  fields.  Also,  the  mathematics  used  in  the  program  are  presented. 

The  initial  computer  program  was  written  for  the  simple  chemical  system 

N2O4  +  Nz1Z12N02  +  N2  .  This  system  was  chosen  so  that  basic  programming 

difficulties  would  not  be  concealed  in  a  maze  of  chemical  reactions  and  also  because  some 
excellent  experimental  data  were  available  from  the  Jet  Propulsion  Laboratory,  which 
could  be  used  for  comparison  purposes.  After  completion,  the  program  was  modified  to 
handle  the  combustion  products  of  hydrogen  and  air  and  to  compare  the  calculations  with 
experimental  data  from  NASA,  Lewis  Research  Center. 

Excellent  agreement  was  obtained  between  the  computed  flow  properties  and  the  experi¬ 
mental  data  for  the  N204  +  2N02  +  N2  system.  Although  some  problems  were 

encountered  in  the  interpretation  of  the  data,  good  agreement  was  obtained  between  the 
computed  data  and  experimental  data  for  the  hydrogen-air  system. 
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SYMBOLS 

speed  of  sound 

characteristic  direction 

specific  heat  at  constant  pressure 

specific  enthalpy  on  common  base 
total  enthalpy 

backward  reaction  rate 

forward  reaction  rate 

concentration  equillbr'um  constant 

Mach  number 

third  body  as  it  applies  to  the  reactions 

molecular  weight 

number  of  species 

number  of  reactions 

pressure 

universal  gas  constant 
temperature 

flow  velocity  In  x  direction 

flow  velocity  in  y  direction 

flow  velocity  tangent  to  streamline 

species  molecular  weight 

variable  cf  rectangular  coordinate  system 

variable  of  rectangular  coordinate  system 

mass  fraction 
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SYMBOLS  (Contd) 

€  constant  in  continuity  equation  (O  for  two-dimensional  and 
1  for  axisymmetric  flow) 

rj  distance  normal  to  streamline  in  natural  coordinate  system 
9  flow  direction 

/x  Mach  angle 

£  distance  along  streamline  in  natural  coordinate  system 

v  stoichiometric  coefficient 

p  density 

<£  fuel-to-air  ratio  over  stoichiometric  fuel-to-air  ratio 

initial  number  of  moles  of  NOg 

SUBSCRIPTS 

i ,  j  refers  to  i  th  and  j  th  species 
f  refers  to  frozen  conditions 

t  total  conditions 

SUPERSCRIPTS 

r'  refers  to  left  side  of  chemical  reaction  equation 
r"  refers  to  right  side  of  rhemical  reaction  equation 
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SECTION  I 
INTRODUCTION 

With  the  Increased  emphasis  being  placed  on  air-breathing  propulsion*  more  accurate 
techniques  must  be  developed  for  analyzing  engine  components.  The  reason  for  this  is 
that*  at  hypersonic  speeds*  small  changes  in  component  performance  can  produce  large 
degradations  in  the  overall  engine  performance.  One  associated  problem  area  is  the  dis¬ 
sociation  of  the  combustion  products  and  the  subsequent  recombination  of  these  products 
in  the  exhaust  nozzle.  With  an  adequate  analytical  technique*  one  may  be  able  to  design 
the  nozzle  to  minimize  the  recombination  losses  while  still  keeping  the  friction  losses 
reasonable  or*  at  the  very  least,  be  able  to  predict  these  losses. 


With  the  introduction  of  high-speed  computers*  the  calculation  of  recombination  losses 
in  the  nozzle  has  received  more  than  its  share  of  attention.  Many  calculations  have  been 
done  on  the  hydrogen-air  system  considering  finite  rate  chemical  reactions.  The  main 
shortcoming  of  these  calculations  was  the  assumption  of  one-dimensional  flow.  The  limi¬ 
tation  of  the  one-dimensional  flow  calculations  is  overcome,  by  this  author*  by  applying 
the  method  of  characteristics  to  a  chemically  reacting  flow  and  developing  an  IBM-7094 
computer  program  to  handle  the  two-dimensional  case  and  the  axiaymmetric  case.  Al¬ 
though  the  calculations  are  more  time-consuming  than  the  restrictive  one-dimensional 
calculations,  a  complete  picture  of  the  entire  nozzle  flow  field  with  chemical  reaction 
is  now  possible. 

In  this  report,  the  problems  encountered  in  developing  the  computer  program  for  the 
two-dimensional  case  and  the  axisymmetric  case  are  discussed.  Also*  the  method-of- 
oharacterl sties  equations  are  derived,  and  the  results  of  the  computer  program  are  pre¬ 
sented. 

Chu  (Reference  1)  presented  the  equations  for  the  method  of  characteristics  for  the 
two-dimensional  and  axisymmetric  cases  including  chemical  reactions.  He  also  3howed 
that  weak  disturbances  in  supersonic  flow  would  propagate  with  the  frozen  speed  of  sound 
even  If  the  flow  was  in  chemical  equilibrium. 

Lick  (Reference  2)  presented  the  equations  of  Reference  1  in  a  form  more  useful  for 
computational  purposes  but  interpreted  the  frozen  speed  of  sound  to  be  that  due  to  only 
the  translational  and  rotational  modes  of  the  molecules. 

Wegener  and  Cole  (Reference  3)  experimentally  proved  that  weak  disturbances  in  super¬ 
sonic  flow  of  a  reacting  gas  do  propagate  with  the  frozen  speed  of  sound.  In  the  theoretical 
calculations  of  the  frozen  speed  of  sound,  Wegener  and  Cole  included  the  contribution  of 
the  vibrational  degrees  of  freedom  In  the  calculated  frozen  speed  of  sound.  The  experiments 
were  performed  very  carefully,  but  the  accuracy  of  the  measured  wave  angles  was  not 
great  enough  to  ascertain  whether  or  not  the  vibration  contribution  to  the  frozen  speed  of 
sound  should  be  included.  Also,  no  mention  was  made  as  to  why  the  vibrational  contribution 
to  the  frozen  speed  of  sound  was  included  even  though  close  agreement  with  the  experi¬ 
mental  results  could  have  been  obtained  if  it  was  not  Included. 

In  Reference  4,  Blend  examined  the  effect  of  the  frequency  of  a  sound  wave  on  the  speed 
with  which  It  will  propagate  through  a  reacting  mixture.  At  low  frequencies,  the  sound  wave 
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propagates  at  the  equilibrium  speed  of  sound  and,  as  the  frequency  is  increased,  the  sound 
wave  travels  at  the  frozen  speed  of  sound,  which  includes  the  vibrational  contribution  tc 
the  speed  of  sound.  As  the  frequency  is  further  increased,  the  wave  travels  at  a  speed 
that  does  not  include  the  vibrational  contribution  to  the  frozen  speed  of  sound. 

From  these  four  works,  it  app  ars  that  there  is  no  general  agreement  as  to  whether  or 
not  the  vibrational  contributions  to  the  frozen  speed  of  sound  should  be  included.  But,  from 
the  data  of  Reference  4,  it  appears  that  one  must  answer  the  question,  what  is  the  frequency 
of  a  Mach  wave,  before  one  will  know  what  to  include  in  the  frozen  speed  of  sound. 

The  present  computer  program  for  the  two-dimensional  and  the  axisymmetric  cases 
does  include  the  vibrational  contribution  to  the  frozen  speed  of  sound. 

SECTION  II 

DERIVATION  OF  METHOD-OF-CHARACTERISTICS  EQUATIONS 

For  the  derivation  of  the  method  of  characteristics  in  two-dimensional  or  axisymmetric 
flow,  the  equations  of  continuity,  momentum,  and  energy  will  be  written  in  the  natural 
coordinate  system  (£,17),  see  Figure  1»  where  £  is  the  distance  along  a  streamline  and 
17  is  the  distance  normal  to  the  streamline. 


Figure  1.  Coordinate  System 
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From  Reference  1  and  2,  one  has  the  following  equations: 
Continuity 


p  d-±  +  V  +  P\i  &  +  P\u  ^  =  o  (1) 

d£  d£  dr)  y 

Momentum 

pv  Jy  +  it  =  o  w 

af  aj 


Energy 


/>V* 


ag 

d£ 


(3) 


(4) 


Rate  of  Species  Change 


V 


dYj 

ae 


where 


Equation  of  State 


P/RT  £ 
i  =  l 


Wi 


djP  .  d/>  dp  +  df  dT  ^  df  aVj 

d£  ’  dP  dZ  dT  d{  dY.  dZ 


(5) 


(6) 


(7) 
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Substituting  Equation  7  into  Equation  1  gives 

P  AX  +  EX  AX  El  Al  tl  1  _L  ill 

d{  p  '  T  df  "  ^  w,  aj 

£Vwi  "l 


+  pv  AH  +  p\k 


where 


•>  =  I  *>i  Yi 
i=i 

AX-.  £  ill  +  £  Y.  iii. 
a£  it',  1  ae  it,  1  ac 

3hj__  dh|  dj  _  dl_ 
di  dl  di  Pi  di 


Combining  Equations  9a,  9b,  and  4  gives 

dl  _  y 

di  "  i  = 


y  _ v  _dv_ 

i=i  cp  di  cp  d£ 


Substituting  Equation  10  into  Equation  8  gives 

p  Ax  +  ±i  el  +  H 


Pm  £ 

n  2. 

?.VWi  i  = 


_£_v  il  + 

Pm 

V  hi  avj  t  ^V1 

dv 

p  d£ 

j 

i  =  1  Cp  di  CpT 

ae 

1  ay, 

,  Wj  ac  + 

Pm 

dd  sin# 

......  +  p  -  u 

377  rve  y 

Eliminating  In  Equation  11  by  use  of  Equation  2,  one  obtains 

di 

J_  dP_  _PM_  dP_  _PM_  "  _hj_  dlj_  V  3p 

V  d{  +  p  a?  +  T  t,  Cp  af  CpT  3£ 

#V  v  I  dYi  dd  ,  „  sin  0 

"  it,  Wj  ae  377  y 

Z  Yj  /Wj 

i=  I 
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Replacing 


ill 

ac 


■  by  Equation  5  gives 
I  <?P  +  PM  dP 


v  p 


Pv 


l 

i=  I 


v  d? 


Cp  TV  1 


CpT  di 


(12) 


n 

-Pwl 


m 


WjV 


y,  +  pm 


dd 

dr) 


One  must  also  utilize  the  following  relations: 

dP  =  —  d£  + 


PVc 


si 


n  0  ~ 


dP 

d-q 


drj 


(13) 


dd  aC  d6  . 

*8  -  +  Triv 


(14) 


Equations  3,  12,  13,  and  14  are  a  system  of  equations  for  the  unknowns  dP/d£, 
dP/dr),d6/d£,dd/d tj  and  have  a  unique  solution  unlesB,  when  the  system  is  solved 
for  one  of  the  unknowns,  the  value  of  determinant  of  the  coefficients  for  the  solution  is 
zero. 

To  reduce  the  problem  to  the  solution  of  ordinary  differential  equations,  we  require 
that  the  equations  be  linearly  dependent  along  some  characteristic  curve.  If  this  is  true, 
then  at  least  one  of  the  variables  is  arbitrary  and  hence  may  be  discontinuous  along  the 
characteristic  curve.  For  one  variable  to  be  discontinuous,  we  solve  for  its  value  and 
set  the  determinant  of  the  coefficient  of  the  numerator  equal  to  zero  since  we  have  assumed 
the  equations  to  be  linearly  dependent  and  the  determinant  of  the  denominator  is  necessarily 
zero. 


Letting  A  =  P\!  /  P  -  l/V- V/CpT  and  setting  up  the  determinant  of  the  coefficients  of 
the  variable,  one  obtains 


dP  dP  dQ  66 
d£  dr/  d£  drj 

A  0  0  PV 

0  I  Pmz  0 

d£  d7)  0  0 


0  =  A  pVz  +  pV d£2 


0  0  d(  di) 


d7?  ,  +  _ l_ 

d£  ~  -/l v" 


V(/jv/p-,/v-  v/cpT;v 
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Since 

°f  =  cprt/[»cp-r] 


then 


d  rj 


+ 


The  rate  of  change  of  77  with  respect  to  £  along  the  characteristic  curve  is  defined  as  the 
tangent  of  the  Mach  angle  and  is  the  slope  of  the  characteristic  curve  or 


Tan 


(15) 


Now  setting  up  the  numerator  determinant  and  putting  its  value  equal  to  zero,  we  can 
determine  our  other  equation.  It  does  not  matter  which  variable  one  solves  for,  since  the 
same  answer  will  be  obtained  for  any  variable.  Choosing  the  determinant  for  dP/drj  and 
letting 


B 


/"i  %  yi  - 


M  Wi 


1=1  vv 


Y.  _  pW( 


sin  6 

y 


gives 


ABO  pV 

0  0  pV*  0 

d£  dP  0  0 

0  6&  d£  dr) 


0 
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Expanding  the  determinant,  one  has 

-  Ad^dP/)Vl  +  df  [b^V*  d*n  -  dd  p*  V3  ]  =  0 


and  dividing  by  pz  V3  6t)  and  changing  signs  gives 

AV 


Since 


then 


dP  +  dd  =  —  dC 

p\!  Cr)  p\! 


/d7?V  _j_ 

\d£/'  AV 


ton 


dP 


X  ± 


dd 


2  vy  ^ 

ton  p^  Pw  tan  ^ 


=  B 


dP 


tan  p 


±  dd  = 


tan 


dC 


^V 


B 


If  dc  is  the  distance  along  the  characteristic  curve,  one  has 

d£  =  dc  cos  p^ 


but 


dP 

pVl  tan  p 


+ 


dd 


sin  pj  dc  B 


therefore 

smM(  =  -jL  --  o,/V 


dP _ 

pM1  tan  pi 


±  dd 


iLdc  f -«_!■!»  -f 
v  L  y  1=1 


Cp  TV 


Y;  +  1 


77? 


1  =  1  w.v 


To  change  from  £  ,  T)  coordinate  system,  reference  is  made  to  a  diagram  of  the  flow 
field,  see  Figure  1.  Sinoe  di^/d^  along  the  characteristic  curve  dc  is  ±  ton^  , 
dy/dx  along  the  characteristic  curve  must  be  tan  (Q±p^).  The  system  of  dif¬ 

ferential  equations  is  now: 


Along  Characteristics 

d  y 

—  =  tan  (£  ±  p  j  ) 
d« 


(17) 
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s 


dP 


pV  ton^,f 


±  d  9 


dc 


[- 


€  S 


n  9 


-I  — ^  Y.  +  Z 


m 


M  C0TV  '  i“|  W.  V 


',]  <18> 


Along  Streamlines 


> 

* 


VdV  +  dP/p  =  0 


(19) 


udy  —  vdx  =  0 


(20) 


VdYj  =  Y|  dC 


(21) 


The  total  enthalpy  Is  assumed  to  be  the  same  everywhere  in  the  flow  field 

Ht  =  h  +  V2/  2 

To  solve  the  equations,  we  must  set  the  equations  up  in  a  finite  difference  form  and  start 
with  an  initial  line  of  data  that  is  not  along  a  characteristic  line.  In  finite  difference  forms, 
the  equations  are  along  the  characteristic 


Ar  r  ,on  [e  ±  ^f] 


Ay 


where  the  bars  over  the  symbols  denote  average  values  between  two  points  (see  Figure  2): 

>  i  n  9 


Ap 


'f 


—  ±  A  £  =  —  Ac 

P V  tan  jl  f  V 


-n  2 


^  €  sin  a  ^ 


h: 


i  =  l  CPTV 


Vi  +X 


i=i  w-,v 


Along  the  streamlines,  we  have 


vAv  +  Ap  =  o 

7 


Ay  =  v 
Ax  u 
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where 


v  ,  J7~7~s 


To  find  the  location  of  the  next  point  on  our  characteristic  from  two  initial  points,  the 
difference  equat  >ns  are  set  up  and  solved  simultaneously  as  follows: 


ton  ( 0  + 


ton  (0  -  /xf)5_2 

where  the  subscripts  3-1  and  3-2  indicate  between  which  two  points  the  average  value  is 
to  be  taken.  Thus 

y*  "  yi  ~  **  ton  (0  -/xf )  +  *i  ton  (0  +  /xf) 

x3  = - - — - - - - — - - 

ton (0  +  iif)  -  ton  (0  -lit) 

'  '  3  -  I  '3-2 


y3  =  yi  +  (*3"  *1  )  ,0n  (0  +  7*f 

Now  letting 

A,  =  (p\*  ton  /xf)5  ) 


A2  =  (pV*  ton/uf)j  2 
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t 


B,  =  sin  LLf 

'  r  '  3- 1 


Ac 


3-1 


=  i xs~ x.)*  +  (*3-  y.)* 


B-  =  sin  u.l 

c  r  t3-* 


ac,.2  =  yu-  «2)‘  +  ( »,-  »y 

finally  one  has 


p  -  p 
rs  ri 


+  0,  -  0, 


=  B. 


P  -  P 

r3  r* 


0,  +  0? 


=  B, 


* 


p3  :  Bi  f  bz  +  9f  -  6Z  +  p,/a,+  p2/a2 
l/A,+  |/a2 

ej+  e,  +  e,-  (v=_5 


One  also  knows  that  the  streamline  passing  through  Point  3  starts  from  some  point  between 
1  and  2.  This  point  can  be  found  using  the  following  relations  along  a  streamline 


Ay 

Ax 


v 


u 


vAv  +  A  P/p  =  o 

Assuming  a  location  for  the  streamline  origin,  one  linearly  interpolates  between  the 
conditions  at  1  and  2  to  find  the  conditions  at  the  origin  of  the  streamline.  Then 
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* 


where  3-4  Indicates  that  the  average  value  is  taken  between  Points  3  and  4  of  Figure  2. 
Since  we  know  #3(  we  can  find  u3  and  v3  and  check  to  see  if 


This  is  continued  until  the  preceding  relation  is  satisfied.  Then  from 

=  t~  °"a  =  A-  +  U-7T 

3-4 

one  can  determine  the  mass  fractions  at  3.  All  conditions  at  3  are  now  known.  The  process 
is  then  continued  for  the  entire  flow  field. 

SECTION  m 

INITIAL  COMPUTER  PROGRAM  FOR  THE  N204  +  N^2N02  +  N2  SYSTEM 
PROBLEMS  ENCOUNTERED 

The  initial  computer  program  was  written  for  the  system  N204+N2X 2N0Z  +  N2  so  that 
errors  in  the  characteristic  equations  would  not  be  concealed  in  a  maze  of  chemical  equations. 
The  initial  computations  were  for  the  two-dimensional  case.  The  nozzle  described  in  Ref¬ 
erence  5  was  used  in  the  calculations. 

To  start  the  method-of-characteristics  calculations,  one  must  know  the  conditions  at  various 
points  on  an  initial  data  line  that  is  not  a  characteristic  line.  For  the  two-dimensional  nozzle 
considered,  we  used  one-dimensional  calculations  with  finite  reactions  rates  to  obtain  this 
data.  For  a  point  close  to  the  throat,  the  conditions  on  the  initial  data  line  were  assumed  to 
be  the  same  as  the  one-dimensional  case.  Because  of  the  low  divergence  angle  of  the  nozzle 
in  Reference  5,  this  is  a  good  approximation. 

Also,  since  finite  difference  techniques  are  used,  average  values  of  the  variables  must  be 
used  in  the  equations  and  one  must  assume  initially  that  the  flow  properties  at  the  point  being 
solved  for  are  the  same  as  at  the  initial  point.  One  can  then  solve  forP,  @and  all  other  vari¬ 
ables  at  'he  new  point  and  use  their  new  values  to  find  the  average  values  for  the  next  cal¬ 
culation.  This  procedure  was  then  repeated  until  consecutively  calculated  pressures  agreed 
within  0.1  percent.  The  only  problem  encountered  in  these  calculations  was  that  the  rates  of 
change  of  the  species  mass  fractions  were  so  rapid  that  finite  differences  could  not  be  used 
to  obtain  the  mass  fractions  at  the  new  point. 

Since  one  has  already  calculated,  at  least  to  a  first  approximation,  the  values  of  x,  y,  P, 
6,  and  V  at  the  new  point  in  the  flow  field,  this  problem  can  be  circumvented.  If  one  connect®! 
the  two  initial  points  by  a  straight  line,  the  point  on  this  line,  from  which  the  streamline 
passing  through  the  new  point  originates,  can  be  found.  The  mass  fractions  can  then  be 
found  by  numerical  Integration  of  Equation  21.  To  perform  this  numerical  integration,  one 
must  know  the  temperature,  density,  and  velocity  along  this  approximate  streamline.  These 
can  be  found  by  first  assuming  that  the  specific  heat  and  molecular  weight  at  the  new  point 
are  the  same  as  at  the  original  point,  or 
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(H=  (S), 

m,  --  m, 


Next  one  assumes  that  T ,  p,  and  V  vary  linearly  along  this  streamline  and  one  can  then 
numerically  integrate  the  equation  for  the  change  in  mass  fractions.  This  must  be  done  for 
every  iteration  during  the  search  for  P3  and  is  the  most  time  consuming  part  of  the  program. 

COMPARISON  OF  EXPERIMENTAL  DATA  WITH  CALCULATED  RESULTS 

Since  the  nozzle  of  Reference  5  had  a  very  small  expansion  angle  (~1.15  degrees) ,  one 
would  expect  the  flow  to  be  nearly  one-dimensional  so  that  the  method-of-characteristics 
calculations  should  be  very  close  to  the  one-dimensional  calculations.  As  can  be  seen  In 
Figures  3  through  6,  this  Is  the  case. 

Figure  3  compares  the  experimental  pressure  distribution  for  the  two-dimensional  flow 
with  that  calculated  one -dimensionally  and  by  the  method  of  characteristics.  The  agreement 
Is  quite  good. 

Figure  4  shows  the  theoretical  temperature  distribution  through  the  nozzle,  although  no 
experimental  temperature  measurements  were  made.  One  can  see  that  even  though  the 
pressure  distribution  Indicated  that  the  flow  is  in  equilibrium,  the  temperature  distribution 
indicates  a  slight  departure  from  equilibrium. 

Figure  5  is  a  comparison  of  the  experimental  and  theoretical  concentration  ratio  of  N00. 
The  experimental  accuracy  of  these  measurements  is  not  as  great  as  for  the  pressure 
measurements  but  the  agreement  is  still  satisfactory. 

Figure  6  shows  the  theoretical  distribution  of  the  mass  fractions  of  NOz. 
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Figure  3.  Nozzle  Wall  Pressure  versus  Distance  Using  Two-Dimensional 
Nozzle  Geometry  from  Reference  5 
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Figure  4.  Temperature  at  Nozzle  Wall  versus  Distance  Using  Two-Dimensional 
Nozzle  Geometry  from  Reference  5 
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SECTION  IV 

MODIFICATION  OF  COMPUTER  PROGRAM  FOR  THE  HYDROGEN-AIR  SYSTEM 

REACTIONS  CONSIDERED  AND  THEIR  RATES 

For  the  hydrogen-air  system,  14  reactions  were  assumed  to  be  occurring  among  the 
9  species  considered.  The  14  reactions  were  assumed  to  be  simultaneously  occurring  in 
the  flow  stream.  These  reactions  and  their  rates  are: 

Reaction  No.  Chemical  Reaction 


1 

h2  +  o2  n 

OH  +  OH 

2 

h  +  o2  r 

OH  +  O 

3 

o  +  h2  r 

OH  +  H 

4 

OH  +  h2  r 

HzO  +  H 

5 

OH  +  OH  r 

h2o  +  0 

6 

H  +  H  +  M  n 

H2  +  m 

7 

H  +  OH  +  Mr 

:h2o  +  M 

8 

H  +  0  +  M  r 

OH  +  M 

9 

o  +  o  +  m  r 

02  +M 

10 

N  +  N  +  M  r 

N2  +  M 

11 

N  +  0  +  M  r 

NO  +  M 

12 

n  +  o2  r 

NO  +  0 

13 

n2  +  o  r 

NO  +  N 

14 

n2  +  o2  r 

2NO 

The  M  refers  to  some  third  body, 
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The  reaction  rates  for  reactions  1  through  8  were  taken  from  Reference  6  and  the  rates 
for  reactions  9  through  14  were  obtained  from  Reference  7.  These  rates  are: 


Rate 


Unit 


=  1  x  1014  exp  (-35229/T*K) 

kf2 
^3 
kf4 
kf5 
kf6 

kf7 
kf8 
kf9 

kno 
kni 

kfl2 


3  x  1014  exp  (-8810/T*K) 

3  x  1014  exp  (-4030/T*K) 

3  x  1014  exp  (~3020/T*K) 

3  x  1014  exp  (-3020/T*K) 
,15 


5  x  10 


1  x  10 


17 


*  1  x  10 


16 


3  x  10 


14 


3  x  10 


14 


14 


=  6  x  10 

=  3.2  x  1012(XK)  1/2  exp  (-3120/T°K) 


cm 


mole-sec. 


cm 


mole  -sec. 


cm 


mole-sec. 


and 


.13 


k^g  =  5  x  10  exp  (-37790/T’K) 


.13 


kfl4  =  2,7  X  10  GXp  (~53842/T#K) 


kbi 


-  kfi/Kc 


where  i  refers  to  the  reaction  number. 


In  the  extension  of  the  program  to  the  hydrogen- air  system,  the  calculations  were  for  the 
axisymmetric  case.  Since  the  consideration  of  this  number  of  species  increased  the  numerical 
integration  time  by  a  large  amount  over  that  required  for  the  N204  system,  a  major  effort 
was  made  to  reduce  this  time  to  a  minimum. 

Because  numerical  integration  requires  the  calculation  of  the  reaction  rates  at  a  large 
number  or  subintervals  In  any  given  interval  and,  comparatively  speaking,  exponential 
functions  require  a  long  computation  time,  an  examination  was  made  to  see  which  reaction 
rates  could  be  written  in  terms  of  other  reaction  rates.  By  minor  changes  in  the  exponential 
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factors  In  the  rates  k^,  k^g,  and  k^,  one  can  write  these  reaction  rates  in  terms  of 
and  k^g.  One  now  has 


r  2  14 

kfl  ?  IX  I014  exp  {-  35240  /  T°K)  =  I  X  I014  - [4- 

fl  L  3  X  10  J 


.  ,*  r  kf  "|5 

- S~7r\  - L“ 7T 

3  X  10  J  L  3  X  10  J 


is 


kf|4  =  2.7  X  10  exp 


(-  53890/ T°K)  =  2.7  X  lo'3  [— — f — j^l 

1  3X  10  J  1  3  X  10  J 


One  other  means  of  reducing  the  integration  time  is  to  combine  as  many  terms  as  possible 
in  the  equations  for  the  rates  of  change  of  each  species.  For  the  9  species  and  14  reactions 
considered,  there  are  88  terms  in  the  equations  for  the  rateB  of  change  of  all  the  species, 
each  term  contains  approximately  3  multiplications.  Of  these  88  terms,  only  40  are  independent 
of  the  other  terms  so  that  by  combining  the  terms  one  can  save  192  multiplications  in  each 
integration  subinterval. 

The  next  problem  was  to  find  the  proper  combination  of  specified  accuracy  on  the  pressure 
and  mass  fractions  that  would  produce  accurate  answers  and  the  shortest  running  times. 
It  was  found  that  when  specifying  5  places  accuracy  on  the  pressure  and  mass  fractions,  the 
program  required  15  to  20  iterations  to  converge  on  the  pressure.  When  the  accuracy  on 
the  pressure  was  decreased  to  3  significant  figures,  only  2  to  3  iterations  were  required 
and  the  answers  were  essentially  the  same  for  both  cases.  Also,  suprisingly  enough,  de¬ 
creasing  the  accuracy  on  the  mass  fractions  increased  the  number  of  iterations  required 
to  converge  on  the  pressure.  Thus  5  places  accuracy  was  retained  on  the  mass  fractions. 

PROBLEMS  ENCOUNTERED  IN  STARTING  THE  SOLUTION 

The  next  step  was  to  compare  the  calculations  with  the  experimental  data  of  Reference  8. 
The  first  attempt  at  obtaining  the  conditions  on  the  Initial  data  line  was  to  follow  the  same 
procedure  as  for  the  NgO+  system.  Upon  proceeding  in  this  manner,  disturbances  were 
noticed  in  the  solution  after  only  a  few  pointB  downstream  of  the  initial  data  line.  These 
disturbances  were  carried  through  the  flow  field  grid  until  they  produced  a  discontinuity 
In  the  center  line  pressure  at  a  distance  of  0.1  foot  downstream  of  the  throat.  This  prompted 
a  closer  look  into  the  flow  in  the  throat  region.  It  was  found  that  Reference  9  provided  a 
very  convenient  analytic  solution  for  the  flow  in  this  region.  To  gain  confidence  in  the 
accuracy  of  the  equations  developed  in  Reference  9,  we  computed  the  location  and  shape 
of  the  sonic  line  for  the  cold  flow  calibration  tests  of  Reference  8.  The  surface  area  of 
revolution  of  the  sonic  line  was  fourd  to  be  0.0548  square  foot  whereas  the  geometric  throat 
area  was  0.0517  square  foot.  Using  the  sonic  line  surface  area  and  the  measured  total 
temperature  and  pressure  of  543*R  and  6892  pounds  per  square  foot,  respectively,  we 
computed  the  theoretical  mass  flow  to  be  8.664  pounds  per  second  whereas  the  measured 
mass  flow  was  8.637  pounds  per  second,  a  difference  on  only  0.13  percent. 

The  transonic  equations  were  used  to  compute  the  coordinates  for  the  Mach  number  of 
1.2  line  and  to  start  the  solution  along  this  line.  Because  of  the  parabolic  shape  of  this  line. 
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there  was  a  sudden  change  in  grid  size  at  the  center  line  downstream  of  the  starting  point. 
This  sudden  change  in  grid  size  produced  erroneous  results  at  the  station  at  which  this 
occurred. 

The  next  step  waste  assume  a  straight  initial  data  line  and  then  to  compute  the  Mach  number 
on  this  line  at  11  equally  spaced  points.  The  initial  data  at  each  of  these  points  was  taken  to 
be  the  same  as  that  obtained  from  the  one-dimensional  analysis  at  the  same  Mach  number  as 
that  at  the  point  of  the  data  line.  This  eliminated  the  problem  of  the  sudden  grid  size  change 
and  produced  smooth  curves. 

COMPARISON  OF  EXPERIMENTAL  DATA  WITH  CALCULATED  RESULTS 


The  comparison  of  the  theoretical  and  experimental  results  for  the  hydrogen-air  system 
are  presented  in  Figures  7  through  15.  The  computed  results  are  compared  with  the  experi¬ 
mental  results  obtained  from  the  nozzle  of  Reference  8.  The  physical  measurements  of  the 
nozzle  are  shown  in  Figure  7.  The  locations  of  the  static  pressure  taps  are  shown  in  the 
throat  regpn  along  with  the  optical  port  locations  at  which  static  temperature  measurements 
were  made. 


Figure  8  shows  the  nozzle  contour  in  the  thro  at  region  and  gives  the  equation  used  to  describe 
this  contour. 

Figure  9  shows  a  portion  of  the  characteristic  net  generated  by  the  program. 


Figure  10  shows  a  comparison  of  the  pressure  distribution  obtained  by  the  method  of  charac¬ 
teristics  with  the  experimental  data  and  with  a  one-dimensional  analysis.  One  of  the  main 
difficulties  in  performing  this  comparison  was  in  obtaining  the  location  of  the  static  taps  as 
is  Indicated  by  the  bar  lines.  The  “kink”  in  the  wall  pressure  distribution  obtained  by  the 
method  of  characteristics  is  due  to  the  discontinuity  in  the  second  derivative  of  the  equations 
used  to  describe  the  wall  shape;  it  will  occur  in  any  similar  case  even  though  the  slope  of 
the  wall  is  continuous.  This  discontinuity  will  produce  a  shock  wave  downstream  of  this  position 
in  the  theoretical  analysis,  but  In  the  physical  situation  this  discontinuity  appears  to  be 
smoothed  over  by  the  boundary  layer  (Reference  10). 

Figure  11  shows  the  computed  temperature  distribution  along  with  the  experimental  points 
and  the  one-dimensional  temperature  distribution.  Since  the  experimental  temperatures  were 
determined  by  the  sodium-D  line  reversal  method,  one  does  not  obtain  a  temperature  at  the 
wall  or  centerline,  but  an  average  temperature  of  the  stream.  Therefore,  the  only  conclusion 
one  can  reach  for  an  axlsymmetric  nozzle  is  that  the  measured  temperature  should  be  some¬ 
where  between  the  calculated  stream  temperature  at  the  wall  and  at  the  centerline. 

Figures  12  through  15  show  comparisons  of  the  mass  fractions  of  H,  OH,  H2 ,  and  H20 
obtained  by  the  method  of  characteristics  with  those  calculated  one-dimensionallv. 
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Figure  7.  Dimensions  of  Nozzle  Used  for  Axisymmetric  Flow  Field  Calculations 
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Figure  10.  Comparison  of  Theoretical  Pressure  Distribution  with 
Experimental  Data 


2- 


5210 


AFAPL-TR-65-20 


(  yo )  ejntojsdujoi 


25 


Figure  11.  Comparison  of  Theoretical  Temperature  Distribution 
Experimental  Data 
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Figure  12.  Variation  in  Mass  Fraction  of  H  for  Axisymmetric  Nozzle 
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Figure  13.  Variation  in  Mass  Fraction  of  OH  for  Axisymmetric  Nozzle 
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Figure  15.  Variation  in  Mass  Fraction  of  H?0  for  Axisymmetrie  Nozzle 
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SECTION  V 
CONCLUSIONS 

Through  the  use  of  a  recently  developed  computer  program,  the  method  of  characteristics 
has  successfully  been  applied  to  the  computation  of  the  flow  field  of  a  chemically  reacting 
gas  in  a  DeLaval  nozzle.  The  accomplishment  enables  one  to  solve  two-dimensional  and 
axisymmetric  reactive  flow  fields. 

Excellent  agreement  was  obtained  between  the  computed  flow  properties  and  the  experi¬ 
mental  data  for  the  NpO  +  N  ZZ  2NO„  +  N,  system.  Good  agreement  waB  obtained  with  the 
experimental  pressure  clata  ?or  the  hydrogen-air  system.  No  direct  comparison  with  the 
temperature  measurements  of  Reference  8  was  possible  since  the  temperatures  measured 
were  average  temperatures  and  not  specifically  the  centerline  or  wall  temperature.  However, 
at  the  second  temperature  measuring  station  (Figure  7),  the  calculated  temperature  profile 
was  fairly  uniform  and  agreed  quite  well  with  the  measured  temperature. 

In  addition  to  analyzing  chemically  reacting  nozzle  flows,  one  can  use  the  program  for 
examining  the  combustion  of  a  premixed  supersonic  hydrogen-air  mixture.  The  program 
also  provides  a  useful  standard  that  can  be  used  to  compare  approximate  methods  such  as 
those  developed  in  Reference  11. 
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